One major task for scientists in his their daily routine is to prepare compelling graphics for technical document, reports or manuscripts for publication. Graphics in form of figures carry the weight of the arguments. They need to be clear, attractive and convincing. according to [@], the difference between good and bad figures always lead into the difference between a highly influential or an obscure paper; a grant/contract won or lost; and a job interview gone well or poorly.

For the last ten years I have squeezed myself into preparing figures for scientific publications and have made throusands of figures. Honestly to say that over this period I have switched from one software to the other in the figure preparation pipeline. I made figures using Microsfoft Excel, OriginPro, SPSS, Matlab, SigmaPlot, matplotlib in python, base R, ggplot2 in R and many others. However, my current preffered tool for making graphics is the ggplot2 package in R. However, looking on the spectrum over the last ten years and the constant switch from one software/tools to the other, I dont expect that I will continue using ggplot2 for the next ten years.

One thing I have learned over the years is that automation should be the main skills for scientists. Automation serve time for data preparation, analysis and producing outputs. In this post I will take you through how to visualize vector field of surface current from drifter observation using the ggplot2 package [@ggplot]. I will further highlight the drawbacks of the default geom of ggplot2 and why it sometimes fail to produce elegant oceanographic plot. Last I will show you how to use altenative geoms from metR package [@metr] to make overcome the challenges inherited in ggplot2 package.

require(metR)
## Loading required package: metR
require(tidyverse)
## Loading required package: tidyverse
## -- Attaching packages --------------------------------------- tidyverse 1.2.1 --
## v ggplot2 3.1.0     v purrr   0.2.5
## v tibble  1.4.2     v dplyr   0.7.8
## v tidyr   0.8.2     v stringr 1.3.1
## v readr   1.3.0     v forcats 0.3.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
require(lubridate)
## Loading required package: lubridate
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
require(oce)
## Loading required package: oce
## Loading required package: gsw
## Loading required package: testthat
## 
## Attaching package: 'testthat'
## The following object is masked from 'package:dplyr':
## 
##     matches
## The following object is masked from 'package:purrr':
## 
##     is_null
## 
## Attaching package: 'oce'
## The following object is masked from 'package:metR':
## 
##     coriolis
require(ocedata)
## Loading required package: ocedata
require(sf)
## Loading required package: sf
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

etopo bathmetry

etopo = sp::read.asciigrid("e:/GIS/ROADMAP/Etopo 1/Tanzania_etopo1/tanz1_-3432.asc")

etopo = etopo %>% as.tibble() %>% rename(lon = 2, lat = 3, bath = 1) %>% filter(bath >= -2400 & bath < 0)


ggplot() +
# metR::geom_contour2(data = etopo, aes(x = lon, y = lat, z = bath),  breaks = seq(-200,0,50), col = "grey70" )+
# metR::geom_text_contour(data = etopo, aes(x = lon, y = lat, z = bath),breaks = seq(-200,0,50), check_overlap = TRUE, rotate = TRUE)+
# metR::geom_contour2(data = etopo, aes(x = lon, y = lat, z = bath),
#                     breaks = seq(0,4000,200))+
  metR::geom_contour_fill(data = etopo, aes(x = lon, y = lat, z = bath),
                          breaks = seq(-4000,0,400), na.fill = TRUE)+
  metR::geom_contour_tanaka(data = etopo, aes(x = lon, y = lat, z = bath),
                            breaks = seq(-4000,0,400))+
  metR::geom_text_contour(data = etopo, aes(x = lon, y = lat, z = bath),
                        breaks = seq(-4000,0,400), check_overlap = TRUE,
                        rotate = TRUE, size = 3)+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"),
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+
  labs(x = NULL, y = NULL)
## Warning: imputing missing values

Dataset

The drifter dataset contain surface current information linked to locations in the physical world. This spatial information help us to understand where the high speed current versus low speed current are found in the ocean. It is helpful to visualize this kind of data in the proper geospatial context i.e to show the data on a realistic map. I have filtered the data to cover the western part of the tropical indian ocean. I prepared and arrange drifter information in data frame—a rectangular collection of variables (columns) and observations (rows). The dataset contains observation of surface current worlwide collected by the Global Drifter Program on all major oceans. Among the variables in the dataset are shown in table @ref(tab:tab1) include:

The drifter observations were randomly distributed within the area as shown in figure @ref(fig:fig1) and requires binning—a process of making equal size grid in the area.

ggplot() +
  geom_point(data = drifter.kimbiji %>% sample_frac(1) , aes(x = lon, y = lat))+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+  
  labs(x = NULL, y = NULL)
The distribution of drifter observation within the area

The distribution of drifter observation within the area

To minimize biasness of sampling, the area was divided into equal size grids showin in figure @ref(fig:fig2).

## convert drifter observation into simple features
drifter.kimbiji.sf = drifter.kimbiji %>% st_as_sf(coords = c("lon", "lat")) %>%
  st_set_crs(4326)

## divide the tropical indian ocean region into grids
drifter.grid = drifter.kimbiji.sf %>% st_make_grid(n = c(50,50))%>%st_sf()

## 
ggplot()+
  geom_sf(data = drifter.grid)+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))
Gridded area to fill drifter observations

Gridded area to fill drifter observations 

## although coordinates are longitude/latitude, st_contains assumes that they are planar
## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant
## over geometries of x
## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for
## longitude/latitude data

Once the area was gridded, then the the mean value of U and V component and the number of observation were calculated in each grid cell.

drifter.kimbiji.sf.se = drifter.kimbiji.sf%>% filter(season=="SE") 

drifter.gridded = drifter.grid %>% 
  mutate(id = 1:n(), contained = lapply(st_contains(st_sf(geometry),drifter.kimbiji.sf.se),identity),
         obs = sapply(contained, length),
         u = sapply(contained, function(x) {median(drifter.kimbiji.sf.se[x,]$u, na.rm = TRUE)}),
         v = sapply(contained, function(x) {median(drifter.kimbiji.sf.se[x,]$v, na.rm = TRUE)}),
         sst = sapply(contained, function(x) {median(drifter.kimbiji.sf.ne[x,]$sst, na.rm = TRUE)})) 
## although coordinates are longitude/latitude, st_contains assumes that they are planar
## Then convert the gridded drifter information into data frame

drifter.gridded = drifter.gridded %>% select(obs, u, v,sst) 

## obtain the centroid coordinates from the grid as table
coordinates = drifter.gridded %>% 
  st_centroid() %>% 
  st_coordinates() %>% 
  as_tibble() %>% 
  rename(x = X, y = Y)
## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant
## over geometries of x
## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for
## longitude/latitude data
## remove the geometry from the simple feature of gridded drifter dataset
st_geometry(drifter.gridded) = NULL

## stitch together the extracted coordinates and drifter information int a single table for SE monsoon season
current.gridded.se = coordinates %>% bind_cols(drifter.gridded) %>% mutate(season = "SE")

drifter.current.gridded = current.gridded.ne %>% bind_rows(current.gridded.se)

After binning, we found that some grids lack the drifter information, therefore, these grids were filled with the U, V values using the Interpolation technique.

## select grids for SE season only
drf.se = drifter.current.gridded %>%select(-sst)%>%filter(season == "SE") %>% na.omit()

## interpolate the U component
u.se = interpBarnes(x = drf.se$x, y = drf.se$y, z = drf.se$u)
dimension = data.frame(lon = u.se$xg, u.se$zg) %>% dim()

## make a U component data table from interpolated matrix
u.tb = data.frame(lon = u.se$xg, u.se$zg) %>% gather(key = "lata", value = "u", 2:dimension[2]) %>% mutate(lat = rep(u.se$yg, each = dimension[1])) %>% select(lon,lat, u) %>% as.tibble()

## interpolate the V component
v.se = interpBarnes(x = drf.se$x, y = drf.se$y, z = drf.se$v)

## make the V component data table from interpolated matrix
v.tb = data.frame(lon = v.se$xg, v.se$zg) %>% gather(key = "lata", value = "v", 2:dimension[2]) %>% mutate(lat = rep(v.se$yg, each = dimension[1])) %>% select(lon,lat, v) %>% as.tibble()

## interpolate the observations
ob.se = interpBarnes(x = drf.se$x, y = drf.se$y, z = drf.se$obs)

## make the V component data table from interpolated matrix
obs.tb = data.frame(lon = ob.se$xg, ob.se$zg) %>% gather(key = "lata", value = "obs", 2:dimension[2]) %>% mutate(lat = rep(ob.se$yg, each = dimension[1])) %>% select(lon,lat, obs) %>% as.tibble()

## stitch now the V component intot the U data table and compute the velocity
uv.se = u.tb %>% bind_cols(v.tb %>% select(v), obs.tb %>% select(obs)) %>% mutate(vel = sqrt(u^2+v^2))

Visualising SE monsoon

ggplot() +
  metR::geom_contour_fill(data = uv.se, aes(x = lon, y = lat, z = obs), na.fill = TRUE, bins = 40) + 
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-7,-3.4), xlim = c(38.5, 42))+
  scale_fill_gradientn(name = "Current",colours = oceColorsJet(120), 
                       # limits = c(0,16),
                       # breaks =seq(4,16,4),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))
ggplot() +
  metR::geom_contour_fill(data = uv.se, aes(x = lon, y = lat, z = vel), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = uv.se, aes(x = lon, y = lat, dx = u, dy = v), 
                    arrow.angle = 20, arrow.type = "open", arrow.length = .5, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_fill_gradientn(name = "Current",
                       colours = oceColorsVelocity(120), 
                       # limits = c(0,1.6), 
                       # breaks =seq(0.2,1.6,.3), 
                       na.value = "white")+
  labs(x = "", y = "")+
 theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_mag(max = 1.8, name = expression(speed~(ms^{-1})))
vector field showing the surface current speed and direction

vector field showing the surface current speed and direction 

ggplot()+
metR::geom_streamline(data = uv.se, 
                        aes(x = lon, y = lat, dx = u, dy = v, 
                            color = sqrt(..dx..^2 + ..dy..^2), 
                            alpha = ..step..),
                      L = 2, res = 2, n = 60, 
                      arrow = NULL, lineend = "round")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_color_viridis_c(guide = "none")+
  scale_size(range = c(0, 1.5), guide = "none") +
  scale_alpha(guide = "none") +
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+  
  labs(x = "", y = "")
Massless current flow during the SE season

Massless current flow during the SE season 

## select grids for SE season only
drf.se.sst = drifter.current.gridded %>%filter(season == "SE") %>% na.omit()

## interpolate the U component
sst.se = interpBarnes(x = drf.se.sst$x, y = drf.se.sst$y, z = drf.se.sst$sst)
dimension = data.frame(lon = sst.se$xg, sst.se$zg) %>% dim()

## make sst into a long format data table from interpolated matrix
sst.se.tb = data.frame(lon = sst.se$xg, sst.se$zg) %>% 
  gather(key = "lata", value = "sst", 2:dimension[2]) %>% 
  mutate(lat = rep(sst.se$yg, each = dimension[1])) %>% 
  select(lon,lat, sst) %>% as.tibble()
ggplot() +
  metR::geom_contour_fill(data = sst.se.tb %>% filter(sst > 27 & sst <=30), aes(x = lon, y = lat, z = sst), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = uv.se, aes(x = lon, y = lat, dx = u, dy = v), 
                    arrow.angle = 30, arrow.type = "open", arrow.length = .5, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.0,-4.2), xlim = c(38.7, 40.5))+
  scale_fill_gradientn(name = "Temperature",colours = oceColors9A(120), 
                       # limits = c(0,23),
                       # breaks =seq(4,20,4),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+
  scale_mag(name = "Speed", max = 1.8, labels = "1.8 m/s")
## Warning: imputing missing values

Visualising NE monsoon

## nelect grids for ne neason only
drf.ne = drifter.current.gridded%>% select(-sst) %>%filter(season == "NE") %>% na.omit()

## interpolate the U component
u.ne = interpBarnes(x = drf.ne$x, y = drf.ne$y, z = drf.ne$u)
dimension = data.frame(lon = u.ne$xg, u.ne$zg) %>% dim()

## make a U component data table from interpolated matrix
u.tb = data.frame(lon = u.ne$xg, u.ne$zg) %>% gather(key = "lata", value = "u", 2:dimension[2]) %>% mutate(lat = rep(u.ne$yg, each = dimension[1])) %>% select(lon,lat, u) %>% as.tibble()

## interpolate the V component
v.ne = interpBarnes(x = drf.ne$x, y = drf.ne$y, z = drf.ne$v)

## make the V component data table from interpolated matrix
v.tb = data.frame(lon = v.ne$xg, v.ne$zg) %>% gather(key = "lata", value = "v", 2:dimension[2]) %>% mutate(lat = rep(v.ne$yg, each = dimension[1])) %>% select(lon,lat, v) %>% as.tibble()

## interpolate the obnervations
ob.ne = interpBarnes(x = drf.ne$x, y = drf.ne$y, z = drf.ne$obs)

## make the V component data table from interpolated matrix
obs.tb = data.frame(lon = ob.ne$xg, ob.ne$zg) %>% gather(key = "lata", value = "obs", 2:dimension[2]) %>% mutate(lat = rep(ob.ne$yg, each = dimension[1])) %>% select(lon,lat, obs) %>% as.tibble()

## stitch now the V component intot the U data table and compute the velocity
uv.ne = u.tb %>% bind_cols(v.tb %>% select(v), obs.tb %>% select(obs)) %>% mutate(vel = sqrt(u^2+v^2))
ggplot() +
  metR::geom_contour_fill(data = uv.ne, aes(x = lon, y = lat, z = obs), na.fill = TRUE, bins = 40) + 
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_fill_gradientn(name = "Current",colours = oceColorsJet(120), 
                       limits = c(0,23),
                       breaks =seq(4,20,4),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))
ggplot() +
  metR::geom_contour_fill(data = uv.ne, aes(x = lon, y = lat, z = vel), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = uv.ne, aes(x = lon, y = lat, dx = u, dy = v), 
                    arrow.angle = 20, arrow.type = "open", arrow.length = .5, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_fill_gradientn(name = "Current",
                       colours = oceColorsVelocity(120), 
                       limits = c(0,1.6),
                       breaks =seq(0.2,1.6,.3),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_mag(max = 1.6, name = expression(speed~(ms^{-1})))
Vector field superimposed on surface current velocity showing the speed and direction

Vector field superimposed on surface current velocity showing the speed and direction 

ggplot()+
metR::geom_streamline(data = uv.ne, 
                        aes(x = lon, y = lat, dx = u, dy = v, 
                            color = sqrt(..dx..^2 + ..dy..^2), 
                            alpha = ..step..),
                      L = 2, res = 2, n = 60, 
                      arrow = NULL, lineend = "round")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_color_viridis_c(guide = "none")+
  scale_size(range = c(0, 1.6), guide = "none") +
  scale_alpha(guide = "none") +
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+  
  labs(x = "", y = "")
Massless current flow during the NE monsoon

Massless current flow during the NE monsoon 

## select grids for NE season only
drf.ne.sst = drifter.current.gridded %>%filter(season == "NE") %>% na.omit()

## interpolate the U component
sst.ne = interpBarnes(x = drf.ne.sst$x, y = drf.ne.sst$y, z = drf.ne.sst$sst)
dimension = data.frame(lon = sst.ne$xg, sst.ne$zg) %>% dim()

## make sst into a long format data table from interpolated matrix
sst.ne.tb = data.frame(lon = sst.ne$xg, sst.ne$zg) %>% 
  gather(key = "lata", value = "sst", 2:dimension[2]) %>% 
  mutate(lat = rep(sst.ne$yg, each = dimension[1])) %>% 
  select(lon,lat, sst) %>% as.tibble()
ggplot() +
  metR::geom_contour_fill(data = sst.ne.tb %>% filter(sst > 26.5), 
                          aes(x = lon, y = lat, z = sst), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = uv.ne, aes(x = lon, y = lat, dx = u, dy = v), 
                    arrow.angle = 30, arrow.type = "open", arrow.length = .5, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.0,-4.2), xlim = c(38.7, 40.5))+
  scale_fill_gradientn(name = "Temperature",colours = oceColors9A(120), 
                       # limits = c(0,23),
                       # breaks =seq(4,20,4),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+
    scale_mag(name = "Speed", max = 1.8, labels = "1.8 m/s")
## Warning: imputing missing values

sst.anomaly = current.gridded.se %>%select(-season) %>% 
  bind_cols(current.gridded.ne %>% select(sst.ne = sst)) %>% 
  mutate(sst.diff = sst.ne - sst) %>% 
  filter(!is.na(sst.diff))

# interpolate the U component
sst.an = interpBarnes(x = sst.anomaly$x, y = sst.anomaly$y, z = sst.anomaly$sst.diff)
dimension = data.frame(lon = sst.an$xg, sst.an$zg) %>% dim()

## make sst into a long format data table from interpolated matrix
sst.an.tb = data.frame(lon = sst.an$xg, sst.an$zg) %>% 
  gather(key = "lata", value = "sst", 2:dimension[2]) %>% 
  mutate(lat = rep(sst.an$yg, each = dimension[1])) %>% 
  select(lon,lat, sst) %>% as.tibble()
ggplot() +
  metR::geom_contour_fill(data = sst.an.tb, 
                          aes(x = lon, y = lat, z = sst), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = uv.ne, aes(x = lon, y = lat, dx = u, dy = v), 
                    arrow.angle = 30, arrow.type = "open", arrow.length = .5, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.0,-4.2), xlim = c(38.7, 40.5))+
  scale_fill_gradientn(name = "Temperature",colours = oceColors9A(120), 
                       # limits = c(0,23),
                       # breaks =seq(4,20,4),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+
    scale_mag(name = "Speed", max = 1.8, labels = "1.8 m/s")

Wind

load("e:/Data Manipulation/xtractomatic/quikscat_wio_2006.RData")


wind.x = wind_x[["data"]]
wind.y = wind_y[["data"]]
date = wind_x$time %>% as.Date()
lon = wind_x$longitude
lat = wind_x$latitude
wind.x.tb = list()
wind.y.tb = list()

for (j in 1:length(date)){
wind.x.tb[[j]] =  data.frame(lon, wind.x[,,j] %>% as.data.frame()) %>%
  gather(key = "lata", value = "wind_x", 2:362) %>%
  mutate(lat = rep(lat, each = 321), time = date[j]) %>%
  filter(lon >=38 & lon <= 42 & lat >=-8 & lat <= -3)

wind.y.tb[[j]] =  data.frame(lon, wind.y[,,j] %>% as.data.frame()) %>%
  gather(key = "lata", value = "wind_y", 2:362) %>%
  mutate(lat = rep(lat, each = 321), time = date[j])%>%
  filter(lon >=38 & lon <= 42 & lat >=-8 & lat <= -3)
}
wind.x.tb = wind.x.tb %>% bind_rows()
wind.y.tb = wind.y.tb %>% bind_rows()

wind.data = wind.x.tb %>% bind_cols(wind.y.tb %>% select(wind_y))
wind.data = wind.x.tb %>% 
  bind_cols(wind.y.tb)  %>%
  mutate(month = month(time),
         season = month,
         season = replace(season,season %in% c(11,12,1,2,3, 4), "NE"),
         # season = replace(season,season %in% c(4,10), "IN"),
         season = replace(season,season %in% c(5,6,7,8,9,10 ), "SE"))%>%
  select(lon, lat, time,season, wind_x, wind_y)
wind.data.season = wind.data %>% 
  group_by(lon,lat, season) %>%
  summarise(wind.x = mean(wind_x, na.rm = TRUE),
            wind.y = mean(wind_y, na.rm = TRUE),
            speed = sqrt(wind.x^2 + wind.y^2))
ggplot() +
  metR::geom_contour_fill(data = wind.data.season%>%filter(season == "SE"), 
                          aes(x = lon, y = lat, z = speed), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = wind.data.season%>%filter(season == "SE"), 
                    aes(x = lon, y = lat, dx = wind.x/10, dy = wind.y/10), 
                    arrow.angle = 30, arrow.type = "open", arrow.length = .75, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_fill_gradientn(name = "Current",
                       colours = oceColorsVelocity(120), 
                       # limits = c(0,1.6), 
                       # breaks =seq(0.2,1.6,.3), 
                       na.value = "white")+
  labs(x = "", y = "")+
 theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_mag( name = expression(speed~(ms^{-1})))
## Warning: imputing missing values

ggplot() +
  metR::geom_contour_fill(data = wind.data.season%>%filter(season == "NE"), 
                          aes(x = lon, y = lat, z = speed), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = wind.data.season%>%filter(season == "NE"), 
                    aes(x = lon, y = lat, dx = wind.x/10, dy = wind.y/10), 
                    arrow.angle = 30, arrow.type = "open", arrow.length = .75, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_fill_gradientn(name = "Current",
                       colours = oceColorsVelocity(120), 
                       # limits = c(0,1.6), 
                       # breaks =seq(0.2,1.6,.3), 
                       na.value = "white")+
  labs(x = "", y = "")+
 theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_mag(name = expression(speed~(ms^{-1})))
## Warning: imputing missing values